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The purpose of this article is numerical verification of the theory of weak turbulence. 
We performed numerical simulation of an ensemble of nonlinearly interacting free grav- 
ity waves (swell) by two different methods: solution of primordial dynamical equations 
describing potential flow of the ideal fluid with a free surface and, solution of the kinetic 
Hasselmann equation, describing the wave ensemble in the framework of the theory of 
weak turbulence. In both cases we observed effects predicted by this theory: frequency 
downshift, angular spreading and formation of Zakharov-Filonenko spectrum ~ cu"'^. 
To achieve quantitative coincidence of the results obtained by different methods, one has 
to supply the Hasselmann kinetic equation by an empirical dissipation term Sdiss modeling 
the coherent effects of white-capping. Using of the standard dissipation terms from op- 
erational wave predicting model ( WAM) leads to significant improvement on short times, 
but not resolve the discrepancy completely, leaving the question about optimal choice of 
Sdiss open. In a long run WAM dissipative terms overestimate dissipation essentially. 



1. Introduction. 

The theory of weak turbulence is designed for statistical description of weakly- nonlinear 
wave ensembles in dispersive media. The main tool of weak turbulence theory is kinetic 
equation for squared wave amplitudes, or a system of such equations. Since the discovery 
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of the kinetic equation for bosons by Nordheim [[T] (see also paper by Peierls [[2]) in the 
context of sohd state physics, this quantum-mechanical tool was applied to wide variety of 
classical problems, including wave turbulence in hydrodynamics, plasmas, liquid helium, 
nonlinear optics, etc. (see monograph by Zakharov, Falkovich and L'vov [31)- Such kinetic 
equations have rich families of exact solutions describing weak-turbulent Kolmogorov 
spectra. Also, kinetic equations for waves have self-similar solutions describing temporal 
or spatial evolution of weak - turbulent spectra. 

However, in our opinion, the most remarkable example of weak turbulence is wind- 
driven sea. The kinetic equation describing statistically the gravity waves on the surface 
of ideal liquid was derived by Hasselmann [H]. Since this time the Hasselmann equation 
is widely used in physical oceanography as foundation for development of wave-prediction 
models such as WAM, SWAN and WAVEWATCH: it is quite special case between other 
applications of the theory of weak turbulence due to the strength of industrial impact. 

In spite of tremendous popularity of the Hasselmann equation, its vahdity and appli- 
cability for description of real wind-driven sea has never been completely proven. It was 
criticized by many respected authors, not only in the context of oceanography. There are 
at least two reasons why the weak-turbulent theory could fail, or at least be incomplete. 

The first reason is presence of the coherent structures. The weak-turbulent theory 
describes only weakly-nonlinear resonant processes. Such processes are spatially extended, 
they provide weak phase and amplitude correlation on the distances significantly exceeding 
the wave length. However, nonlinearity also causes another phenomena, much stronger 
localized in space. Such phenomena - solitons, quasi-solitons and wave collapses are 
strongly nonlinear and cannot be described by the kinetic equations. Meanwhile, they 
could compete with weakly-nonlinear resonant processes and make comparable or even 
dominating contribution in the energy, momentum and wave-action balance. For gravity 
waves on the fluid surface the most important coherent structures are white-cappings (or 
wave-breakings), responsible for essential dissipation of wave energy. 

The second reason limiting the applicability of the weak-turbulent theory is finite size 
of any real physical system. The kinetic equations are derived only for infinite media, 
where the wave vector runs continuous (i-dimensional Fourier space. Situation is different 
for the wave systems with boundaries, e.g. boxes with periodical or reflective boundary 
conditions. The Fourier space of such systems is infinite lattice of discrete eigen-modes. 
If the spacing of the lattice is not small enough, or the level of Fourier modes is not big 
enough, the whole physics of nonlinear interaction becomes completely different from the 
continuous case. We shall call effects caused by a finite size of a system "mesoscopic 
effects". These effects could be important in nature and they certainly should be taken 
in to account by anyone who performs numerical simulation of wave turbulence. 

For these two reasons verification of the weak turbulent theory is an urgent problem, 
important for the whole physics of nonlinear waves. The verification can be done by direct 
numerical simulation of the primordial djTiamical equations describing wave turbulence 
in nonlinear medium. 

So far, the numerical experimentalists tried to check some predictions of the weak- 
turbulent theory, such as weak-turbulent Kolmogorov spectra. For the gravity wave tur- 
bulence the most important is Zakharov- Filonenko spectrum F^^ ~ uj~^ [ At the 
moment, this spectrum was observed in numerous numerical experiments [[5]-[[2U]. 
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The attempts of verification of weak turbulent theory through numerical simulation of 
primordial dynamical equations has been started with numerical simulation of 2D optical 
turbulence [i21j, which demonstrated, in particular, co-existence of weak-turbulent and 
coherent events. 

Numerical simulation of 2D-turbulence of capillary waves was done in [ [6] , [ [7] , and 
[H]. The main results of the simulation consisted in observation of classical regime of 
weak turbulence with spectrum F^^ ~ a;~^^/^, and discovery of non-classical regime of 
"frozen turbulence" , characterized by absence of energy transfer from low to high wave- 
numbers. The classical regime of turbulence was observed on the grid of 256 x 256 points 
at relatively high levels of excitation, while the "frozen" regime was realized at lower 
levels of excitation, or more coarse grids. The effect of "frozen" turbulence is explained 
by mesoscopic effects: sparsity of both exact and approximate resonances. The classical 
regime of turbulence becomes possible due to broadening of resonances by nonlinearity. 
The classical and the frozen regime could coexist. We call this situation "mesoscopic 
turbulence" . 

In fact, the "frozen" turbulence is close to A'ylM regime, when the dynamics of turbu- 
lence is close to the behavior of integrable system [[8]. 

Simulation of the surface gravity waves turbulence for the first time was done simultane- 
ously by Tanaka [[H] and Onorato at al. [HO]- Due to limitations of calculation performance 
of that time computers simulations were limited to dynamical equations and quite short 
simulation times. After that sime many articles developed these pioneering results [[TT]-[ 
[T7] . The nonlinear effects for gravity waves are weaker than for cappilary waves, as a 
result the influence of the mesoscopic effects is stronger. It makes the simulations much 
more difficult. In the experiments mentioned above, the grid was fine enough to resolve 
the spectral tails and observe the weak turbulent Kolmogorov asymptotics Ik ~ k~^, but 
it was too coarse to avoid mesoscopic effects in the area of spectral peak. The finest grid 
2048 X 4096 was used by Tanaka, however in his experiments the spectral peak was posed 
at k = 68, while the level of nonlinearity was quite small (typical steepness n ~ 0.07). 
More over, his experiment was very short in terms of characteristic time of the waves (only 
several tens periods of the leading wave). If he would continue his calculations he would 
observe formation of the weak turbulent Kolmogorov tail as in the work by Onorato at 
al. [[ID]. However his spectral peak was seriously damaged by mesoscopic effects. 

In this article we present the results of simulation of the surface gravity waves turbu- 
lence, namely modelling of swell propagation. All previous experimentalists since Tanaka [ 
in] and Onorato at al. [TU] performed only one type of experiments — solution of the pri- 
mordial dynamical equations. We made two experiments simultaneously. We solved not 
only dynamical equations but also Hasselmann kinetic equation and compared the results. 
We think that our results can be considered as the first attempt of direct verification of 
Hasselmann kinetic equation. 

In the dynamic experiment we used less fine grid than Tanaka (512 x 4096) but we posed 
the spectral peak at k = 300, and our initial steepness was much higher than in the Tanaka 
case (/X ~ 0.15). Due to much higher k and stronger nonlinearity we managed basically to 
suppress the mesoscopic effects. The bulk of energy containing modes satisfies Rayleigh 
distribution, according to the weak-turbulent scenario. The distribution has a heavy tail 
of abnormally intensive harmonics ("oligarchs" according to terminology, introduced in 
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the article [[15]), but they contain no more than 5% of total energy. 

One important point should be mentioned. In our experiments we observed not only 
weak turbulence, but also additional nonlinear dissipation of the wave energy, which could 
be identified as the dissipation due to white-capping. We observed fast broadening of the 
spectra due to multiple harmonics generation. Formation of the broad spectrum with 
strong small scale tails can be explained and formation of quite sharp structures on waves 
crests. Further development of these tails is suppressed by an artificial dissipation used in 
our dynamical experiment. One can say that in such a way we roughly simulated white- 
capping phenomenon (continuous in time dissipation of energy due to multiple acts of 
wave braking on the very edge of wave crest, arresting formation of derivative singularity 
on the sharp wave crest). Of course we cannot simulate directly wave braking phenomenon 
but one can say that our method provide may be rough but reasonable simulation of this 
effect. 

To reach an agreement with dynamic experiments, we had to add to the kinetic equation 
a phenomenological dissipation term Sdiss- In this article we examined dissipation terms 
used in the operational wave-prediction models WAM Cycle 3 and WAM cycle 4 (hereafter 
referenced as WAM3 and WAM4 correspondingly). Both of these terms overestimate 
nonlinear dissipation significantly. Term given in WAM3 gives acceptable results on short 
periods of time (time less than IOOOTq, where Tq is the time period of the leading wave 
of the initial condition). But at the end of simulation {t = 3378To) an error in wave 
action become close to 30%. For WAM4 the situation is even worse. If the characteristic 
wave length of initial conditions is equal to 22m, then 3378To is a little bit more than 
only 3 hours of wave development. Even primitive viscous dissipative term without any 
simulation of strongly nonlinear events gives us better results. It means that the question 
about a reasonable formulae for Sdiss is open for now. 

2. Deterministic and statistic models. 

In the "dynamical" part of our experiment surface of the fluid was described by two 
functions of horizontal variables x,y and time t: surface elevation r]{x,y,t) and velocity 
potential on the surface ip{x,y,t). They satisfy the canonical equations [ 1^ 



Hamiltonian H is presented by the first three terms in expansion on powers of nonlinearity 



drj 5H dip 
dt^~6^' 'dt 



6H 




(1) 



H 



Ho + H^ + H2 + ..., 




(2) 



Here k is the linear integral operator k = V—V^, defined in Fourier space as 




(3) 
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Using Hamiltonian (j2]) and equations ([T]) one can get the dynamical equations [[6]: 

f] = kip — iy {r]^/ ip)) — k[7]kip\ + 
+k{r)k[r]kilj]) + lV'^[r)'^kil:] + 

lk[rj^V^iP]+F-^[^kVk], (4) 

-[kiP]k[rjkiP] - [r]kij]V^ij + F-^[-fkiJk]. 

Here corresponds to inverse Fourier transform. We introduced artificial dissipative 
terms ^"^[7^7/^] and F~^['jk4'k], corresponding to pseudo- viscous high frequency damping 
following recent work [ [19] . 

The model ([I])-(11D was used in the numerical experiments [E] - [[S], [[12], 

Introduction of the complex normal variables 



k 



where Uk = y/gk, transforms equations ([T]) into 

^ = -i— . (6) 

dt Sat ^ ^ 

k 

To proceed with statistical description of the wave ensemble, first, one should perform 
the canonical transformation —>■ 6^, which excludes the cubical terms in the Hamilto- 
nian. The details of this transformation can be found in the paper by Zakharov (1999) [ 
After the transformation the Hamiltonian takes the form 



H = J u;^,b^,btdk + -J TrkAhWkXK X (7) 
x%+fc^_fc^_fc3dfcidA?2d4. 

where T is a homogeneous function of the third order: 

T{ek, eki, eh, eh) = e^T{k, h, h, h). (8) 

Connection between and 6^ together with explicit expression for T^j:^^^%.^ can be found, 
for example, in [i25j. 

Let us introduce the pair correlation function 

<a^ai >=gN^^5{k-k'), (9) 

where is the spectral density of the wave function. This definition of the wave action 
is common in oceanography. 

We also introduce the correlation function for transformed normal variables 

<h^hl>=gn-^5{k-k') (10) 
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Functions and Nj^ can be expressed through each other in terms of cumbersome power 
series of expansion on /i. Here /i is the characteristic steepness defined as follows 

/^ = |^v^, (11) 

where E is the wave energy and is the wave action. Following this definition for the 
Stokes wave of small amplitude 

r] = a cos{kx), 
11 ~ ak. 

On deep water their relative difference between and N^: is of the order of /i^ and can 
be neglected (in most cases experimental results shows ~ 0.1). 
Spectrum satisfies Hasselmann (kinetic) equation [H] 

— ^ - Sni[n] + Sdiss + 27fcng, 



dt 



\ (12) 



• ^ki^k2'^ks~^ 



x5 [k + ki — k2 — k^) d/cid/c2d/c3. 



Here Sdiss is an empiric dissipative term, corresponding to white-capping. 
Stationary conservative kinetic equation 

Snl = (13) 

has the rich family of Kolmogorov-type [ [26] exact solutions. Among them is Zakharov- 
Filonenko spectrum [ '5] for the direct cascade of energy 

rik ~ ^, (14) 

and Zakharov- Zaslavsky [ES] spectra for the inverse cascade of wave action 

1 , , 



3. Deterministic Numerical Experiment. 

3.1. Problem Setup 

The dynamical equations (j4]) have been solved in the real-space domain 2tt x 27r on 
the grid 512 x 4096 with the gravity acceleration set to g = 1. The solution has been 
performed by the spectral code, developed in [[22] and previously used in [[23], [[12], [[T3].[ 
[15] . We have to stress out that in the current computations the resolution in F-direction 
(long axis) is better than the resolution in X-direction by the factor of 8. 

This approach is reasonable if the swell is essentially anisotropic, almost one-dimensional. 
This assumption will be validated by the proper choice of the initial data for computation. 
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As the initial condition, we used the Gaussian-shaped distribution in Fourier space (see 



The initial phases of all harmonics were random. The average steepness of this initial 
condition was ~ 0.15 (defined in accordance with Eq llll) . 

To realize similar experiment in the laboratory wave tank, one has to to generate the 
waves with wave-length 300 times less than the length of the tank. The width of the tank 
would not be less than 1/8 of its length. The minimal wave length of the gravitational wave 
in absence of capillary effects can be estimated as Xmin — ^cm. The leading wavelength 
should be higher by the order of magnitude A ~ 30cm. 

In such big tank of 200 x 25 meters experimentators can observe the evolution of the 
swell until approximately 700To - still less than in our experiments. In the tanks of smaller 
size, the effects of discreetness the Fourier space will be dominating, and experimentalists 
will observe either "frozen" , or "mesoscopic" wave turbulence, qualitatively different from 
the wave turbulence in the ocean. 

To stabilize high-frequency numerical instability, the damping function has been chosen 

as 



kd = 1024,7 = 5.65 x 10"^ 

The simulation was performed until t = 1225, which is equivalent to 3378 Tq, where 
To = 211 1 ^k^ is the period of the wave, corresponding to the maximum of the initial 
spectral distribution. 

3.2. Zakharov-Filonenko spectra 

Like in the previous papers [ [10] , [ [12] , [ [13] and [ [15] , we observed fast formation of the 
spectral tail, described by Zakharov-Filonenko law for the direct cascade ~ k~^ [ [5] 
(see Figl2]). In the agreement with [[I5], the spectral maximum slowly down-shifts to the 
large scales region, which corresponds to the inverse cascade [I27j.[[28]. 

Also, the direct measurement of energy spectrum has been performed during the final 
stage of the simulation, when the spectral down shift was slow enough. This experiment 
can be interpreted as the ocean buoy record - the time series of the surface elevations has 
been recorded at one point of the surface during T^^^oy — 300To. The Fourier transform of 
the autocorrelation function 
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Figure 1. Initial distribution of |a^| 



on /c-plane. 



allows to detect the direct cascade spectrum tail proportional to u'"^ (see FigJSj), well 
known from experimental observations [ [29] . [ [30] . [ |3T] . 

In this paper we had no intention to improve results of the previous papers. Inertial 
interval of the angle-averaged spectra (tu-spectrum also is angle-averaged because it de- 
pends only on frequency u ~ \/k) is limited due to anisotropy of the integration domain. 
In this case we have less than one decade interval and we cannot find an exponent of the 
Kolmogorov tail surely, however this work was already done in several previous articles 
(see for instance [[13]). Here we limit our observations by qualitative correspondance. 

3.3. Is the weak-turbulent scenario realized? 

Presence of Kolmogorov asymptotic in spectral tails, however, is not enough to validate 
applicability of the weak-turbulent scenario for description of wave ensemble. We have also 
be sure that statistical properties of the ensemble correspond to weak-turbulent theory 
assumptions. 

One should stress out that at the beginning of our experiments |arp is a smooth function 
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Figure 2. Angle-averaged spectrum =< |a^p > in a double logarithmic scale. The tail 
of distribution fits to Zakharov-Filonenko spectrum. 



of k. Only the phases of the individual waves are random. As shows numerical simulation, 
the initial condition f|T6l) (see FigH]) does not preserve its smoothness - it becomes rough 
virtually immediately (see FigHj). The picture of this roughness is remarkably preserved 
in many details, even as the spectrum down-shifts as a whole. This roughness does 
not contradict the weak-turbulent theory. According to this theory, the wave ensemble 
is almost Gaussian, and both real and imaginary parts of each separate harmonics are 
not-correlated. However, also according to the weak-turbulent theory, the spectra must 
become smooth after averaging over long enough time of more than periods. Earlier 
we observed such restoration of the smoothness in the numerical experiments of the MMT 
model (see [Ill|,[ll8], [|19] and [ED]). However, in the experiments discussed in the article, 
the roughness still persists and the averaging does not suppresses it completely. It can be 
explained by sparsity of the resonances. 

Resonant conditions are defined by the system of equations: 

k + ki = k2 + k3, ^ ^ 

These resonant conditions define five- dimensional hyper-surface in six-dimensional space 
k, fci, ^2- In a finite system, ([T9l) turns into Diophantine equation. Some solutions of this 
equation are known [|32], [[17]. In reality, however, the energy transport is realized not by 
exact, but approximate resonances, posed in a layer near the resonant surface and defined 
by 

\^^k + t^fci — i^k2 — ^k+ki-k2 I ^ r, (20) 
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Figure 3. Energy spectrum in a double logarithmic scale. The tail of distribution fits to 
asymptotic cu"^. 



where F is a characteristic inverse time of nonlinear interaction. 

In the finite systems k, ki, k2 take values on the nodes of the discrete grid. The weak 
turbulent approach is valid, if the density of discrete approximate resonances inside the 
layer fl20|) is high enough. In our case the lattice constant is Ak = 1, and typical relative 
deviation from the resonance surface 

Auj uj'i, . , ujI 1 ^ _ , ^2]^^ 



2 X IQ- 



~ —Ak = — ^ 

UJ UJ UJ 600 

Inverse time of the interaction F can be estimated from our numerical experiments: wave 
amphtudes change essentially during 30 periods, and one can assume: T/u ^ 10^^ ^ — . 
It means that the condition for the applicability of weak turbulent theory is typically 
satisfied, but the reserve for their validity is rather modest. As a result, some particular 
harmonics, posed in certain "privileged" point of fc-plane could form a network of almost 
resonant quadruplets and realize significant part of energy transport. Amplitudes of these 
harmonics exceed the average level essentially. This effect was described in the article [ 
[15], where such "selected few" harmonics were called "oligarchs". If oligarchs realize most 
part of the energy flux, the turbulence is mesoscopic, not weak. 

3.4. Statistics of the harmonics 

According to the weak-turbulent scenario, statistics of the aj:{t) in any given k should 
be close to Gaussian. It presumes that the PDF for the squared amplitudes is 

Pila^f) ^ le-l'^.H^/^, (22) 
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Figure 4. Surface |agp before (left) and after (right) low-pass filter at the moment of time 
t ^ 67To. 



where D =< |a^p > is the mean square amplitude. 



To check the equation fl22]) . we need to find the way for calculation of D{k). If the 
ensemble is stationary in time, D{k) could be found for any given k by the time averaging. 
In our case, the process is non-stationary, and we have a problem with determination of 
D{k). 

To resolve this problem, we used low-pass filtering instead of time averaging. The 
low-pass filter was taken in the form 



f{n) 
A = 



= e 
0.25 



■(|n|/A)3 

■ iV^/2, A^^ = 4096. 



(23) 



where n is integer vector running i^-space. 

This choice of the low-pass filter preserves the values of total energy, wave action and the 
total momentum within three percent accuracy. The shape of low-pass filtered function 
is presented on FigjH 

Now it is possible to average the PDF over different areas in /c-space. The results for 
two different moments of time t ~ 67.1To and t ~ 925To are presented in Fig 15] and FigJHl 
The thin line gives PDF after averaging over dissipation region harmonics, while bold 
line presents averaging over the non-dissipative area \k\ < kd = 1024. One can see that 
statistics in the last case is quite close to the Gaussian, while in the dissipation region 
it deviates from Gaussian distribution. However, deviation from the Guassianity in the 
dissipation region does not create any problems, since the "dissipative" harmonics do not 
contain any essential amount of the total energy, wave action and momentum. 

One should remember, that the bold lines in the Fig|5] and Figl6] are the results of 
averaging over a million of harmonics. Among them there is a population of "selected 
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Figure 5. Probability distribution function for relative squared amplitudes |ajtp/ < 
lofcp >. t Q7Tq. Thin and bold lines - averaging over dissipation and dissipation- 
free regions of i?-space correspondingly. 



few", or "oligarchs", with the amplitudes exceeding the average value by the factor of 
more than ten times. The oligarchs exist because our grid is still not fine enough. The 
contribution of such oligarchs in our case to the total wave action does not exceed 4%. 
Fifteen leading oligarchs at some moment of time are presented in the Appendix |Al 
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Figure 6. Probability distribution function for relative squared amplitudes la^p/ < 
lofep >. i ~ 9253]) ■ Thin and bold lines - averaging over dissipation and dissipation- 

— * 

free regions of X-space correspondingly. 
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3.5. Two-stage evolution of the swell 

Fig. ITlfTU] demonstrate time evolution of main characteristics of the wave field: wave 
action, energy, characteristic slope and mean frequency. 




Figure 7. Total wave action as a function of time. Solid line — dynamical equations, 
dashed-dotted line — kinetic equation with artificial viscosity, dashed line — kinetic 
equation with damping term, dotted line — kinetic equation with WAM4 damping 

term 



Figl9] should be specially commented. Here and further we define the characteristic 
slope as defined in Eq. [TTJ According to this definition of steepness for the classical 
Pierson-Moscowitz spectrum fi = 0.095. Our initial steepness fi ~ 0.15 exceeds this value 
essentially. 

Observed evolution of the spectrum can be conventionally separated in two phases. On 
the first stage we observe fast drop of wave action, slope and especially energy. Then 
the drop is stabilized, and we observe slow down-shift of mean frequency together with 
angular spreading. Level lines of smoothed spectra in the first and in the last moments 
of time are shown in Fig JTTlfT^ 

Presence of two stages can be understood through the study of the Probability Distri- 
bution Functions (PDFs) of the elevations of the surface. In all figures we compare the 
distribution in experimental results with Gaussian distribution and Tayfun distribution [ 
l33] . The later case is just a first correction to Gaussian distribution due to small nonlin- 
earity. We used explicit form of Tayfun distribution following [[IB]. In the initial moment 
of time the PDF is Gaussian (see Fig JT3|) . No nonlinear interaction is involved, so Tayfun 
distribution does not fit at all. However, very soon intensive super-Gaussian tails appear 
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Figure 8. Total wave energy as a function of time. Solid line — dynamical equations, 
dashed-dotted line — kinetic equation with artificial viscosity, dashed line — kinetic 
equation with WAM3 damping term, dotted line — kinetic equation with WAM4 damping 
term 



(see Fig JT^ . They are well described by Tayfun distribution. Then tails decrease slowly 
(Fig JTSi) . and in the last moment of simulation, when characteristics of the sea are close 
to Peirson-Moscowitz, the statistics is close to Gaussian again (Fig JT6i) . Moderate tails do 
exist and, in a good agreement with Tayfun correction, troughs are more probable than 
crests which in turn can be much larger in absolute value. PDF for rjy — longitudinal 
gradients in the first moments of time is Gaussian (Fig JT7|) . Then in a very short period 
of time strong non-Gaussian tails appear and reach their maximum at t ~ 14To (Fig lTSi) . 
Here Tq = 27r/ ^/ko — period of initial leading wave. Since this moment the non-Gaussian 
tails decrease. In the last moment of simulation they are essentially reduced(Fig IT9l) . 

Fast growth of non-Gaussian tails can be explained by formation of coherent harmonics. 
Indeed, 14To ^ 2ti /{ujq^) is a characteristic time of nonlinear processes due to quadratic 
nonlinearity. Note that the fourth harmonic in our system is strongly decaying, hence we 
cannot see real white caps. 

Figures [2U1I22] present PDFs for gradients in the orthogonal direction. 
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Figure 9. Average wave slope as a function of time. Solid line — dynamical equations, 
dashed-dotted line — kinetic equation with artificial viscosity, dashed line — kinetic 
equation with WAM3 damping term, dotted fine — kinetic equation with WAM4 damping 
term 
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Figure 10. Mean wave frequency as a function of time. Solid line — dynamical equations, 
dashed-dotted line — kinetic equation with artificial viscosity, dashed line — kinetic 
equation with WAM3 damping term, dotted line — kinetic equation with WAM4 damping 
term 
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Figure 12. Final spectrum |a^p. t = 337870. 
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Figure 13. PDF for the surface elevation rj at the initial moment of time, t — 0. 




Figure 14. PDF for the surface elevation rj at the moment of maximum surface roughness. 
t ~ UTq. 
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Figure 15. PDF for the surface elevation rj at some middle moment of time, t ~ 67.1To. 




Figure 16. PDF for the surface elevation t] at the final moment of time, t = 3378To. 
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Figure 18. PDF for {Vri)y at the moment of maximum surface roughness, t ~ 14To. 
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Figure 20. PDF for (V?7)a; at the initial moment of time, t = 0. 
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Figure 22. PDF for {'V'r])x at the final moment of time, t = 3378To. 
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Figures [23)12^ present snapshots of the surface in the initial and final moments of 
simulation. FigI2S] is a snapshot of the surface in the moment of maximal roughness 
T = 4.94 ~ 14To. 
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Figure 23. Surface elevation at the initial moment of time, t = 0. 



4. Statistical numerical experiment 

4.1. Numerical model for Hasselmann Equation 

Numerical integration of kinetic equation for gravity waves on deep water (Hasselmann 
equation) was the subject of considerable efforts for last three decades. The "ultimate 
goal" of the effort - creation of the operational wave model for wave forecast based 
on direct solution of the Hasselmann equation - happened to be an extremely difficult 
computational problem due to mathematical complexity of the Sni term, which requires 
calculation of a three-dimensional integral at every moment of time. 

Historically, numerical methods of integration of kinetic equation for gravity waves exist 
in two "flavors" . 

The first one is associated with works [El], [ES], [EH], [EZj, [EB] and [EH], and is based 
on transformation of 6-fold into 3-fold integrals using 5-functions. Such transformation 
leads to appearance of integrable singularities, which creates additional difficulties in 
calculations of the Sni term. 

The second type of models has been developed in works of [HQ] and [E], [HI] and is 
currently known as Resio- Tracy model. It uses direct calculation of resonant quadruplet 
contribution into Sni integral, based on the following property: given two fixed vectors 
k, ki, another two ^2,^3 are uniquely defined by the point "moving" along the resonant 
curve - locus. 
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Figure 24. Surface elevation at the final moment of time, t = 3378To. 



Numerical simulation in the current work was performed with the help of modified 
version of the second type algorithm. Calculations were made on the grid 71 x 36 points in 
the frequency-angle domain [u, 6] with exponential distribution of points in the frequency 
domain and uniform distribution of points in the angle direction. 

To date, Resio- Tracy model suffered rigorous testing and is currently used with high 
degree of trustworthiness: it was tested with respect to motion integrals conservation in 
the "clean" tests, wave action conservation in wave spectrum down-shift, realization of 
self - similar solution in "pure swell" and "wind forced" regimes (see [ HI] , [ US] , [ SS] ) • 

Description of scaling procedure from dynamical equations to Hasselman kinetic equa- 
tion variables is presented in Appendix [Bl 

4.2. Statistical model setup 

The numerical model used for solution of the Hasselmann equation has been supplied 
with the damping term in three different forms: 

1. Pseudo- viscous high frequency damping (fTTl) used in dynamical equations; 

2. WAM3 viscous term; 

3. WAM4 viscous term; 

Two last viscous terms are the "white-capping" terms, describing energy dissipation by 
surface waves due to white-capping, as used in WAM wave forecasting models (in SWAN 
model the used term has different but close tunable parameters), see [H6|: 

(I-)' (24, 

where k and u are wave number and frequency, tilde denotes mean value; Cds, S and p 
are tunable coefficients; S = k\fH is the overall steepness; SpM = (3.02 x 10~^)^/^ is 
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Figure 25. Surface elevation at the moment of maximum surface roughness, t ~ 14To. 
Gradients are more conspicuous. 



the value of S for the Pierson-Moscowitz spectrum (note that the characteristic steepness 
yU = v^'S'). It is worth to note that according to [|5T] theoretical value of steepness for 
Pierson-Moscovitz spectrum is SpM — (4.57 x 10~^)^'^. It gives us /i ~ 0.095. 
Values of tunable coefficients for WAM3 case are: 

Cds = 2.36 X 10-^ 6 = 0, p = A (25) 

and for WAM4 case are: 

Cds = 4.10 X 10-^ 5 = 0.5, p = 4 (26) 

In all three cases we used as initial data smoothed (filtered) spectra (see FigH]) obtained 
in the dynamical run at the time T^, = 3.65mm = 24.3 ~ 67.1To, which can be considered 
as a moment of the end of the first "fast" stage of spectral evolution. 

The natural question stemming in this point, is why calculation of the dynamical and 
Hasselmann model cannot be started from the initial conditions (|T6l) simultaneously? 

There are good reasons for that: 

As it was mentioned before, the time evolution of the initial conditions fll6p in presence 
of the damping term can be separated in two stages: relatively fast total energy drop 
in the beginning of the evolution and succeeding relatively slow total energy decrease as 
a function of time, see FiglHl We explain this phenomenon by existence of the effective 
channel of the energy dissipation due to strong nonlinear effects, which can be associated 
with the white-capping as we mentioned in the introduction. 

We have started with relatively steep waves fi ~ 0.167. As we see, at that steepness 
white-capping is the leading effect. This fact is confirmed by numerous field and labora- 
tory experiments. From the mathematical view-point, the white-capping is formation of 
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coherent structures - strongly correlated multiple harmonics. The spectral peak is posed 
in our experiments initially at A; ~ 300, while the edge of the damping area ~ 1024. 
Hence, only the second and the third harmonic can be developed, while higher harmonics 
are suppressed by strong dissipation. Anyway, even formation of the second and the third 
harmonic is enough to create intensive non-Gaussian tail of the PDF for longitudinal 
gradients. This process is very fast. In the moment of time T = 14To we see fully de- 
veloped tails. Relatively sharp gradients mimic formation of white caps. Certainly, the 
"pure" Hasselmann equation is not apphcable on this early stage of spectral evolution, 
when energy intensively dissipates. 

As steepness decreases and spectral maximum of the swell down-shifts, the efficiency of 
such mechanism of energy absorption becomes less important. When the steepness value 
drops down to /i ~ 0.1, at approximately T ~ 280To, the white-capping is negligibly 
small. Therefore, we decided to start comparison between deterministic and statistical 
modeling in some intermediate moment of time T ~ GT.lTg. 

5. Comparison of deterministic and statistical experiments. 

5.1. Statistical experiment with pseudo-viscous damping term. 

The first series of statistical experiments has been performed with pseudo- viscous damp- 
ing term (ITTIl . 

FiglT] - [TU] show total wave action, total energy, mean wave slope and mean wave 
frequency as the functions of time. 

FigjST] shows the time evolution of angle-averaged wave action spectra as the functions 
of frequency for dynamical and Hasselmann equations. We see similar down-shift of the 
spectral maximum both in dynamic and Hasselmann equations. The correspondence of 
the spectral maxima amplitudes is not good at all. 

It is quite obvious that the infiuence of the artificial viscosity is not strong enough to 
reach the correspondence of two models. 

5.2. Statistical experiments with WAM3 damping term 

The second series of statistical experiments has been done for the choice of WAM3 
damping term. 

The temporal behavior of total wave action, energy and average wave slope (see FigJT] 
- [TOl) for WAM3 damping term is in better correspondence with dynamical model, than 
in the case of artificial viscosity term. For initial 50 mm duration of the experiment we 
observe decent correspondence between djTiamical and Hasselmann equations. For longer 
time the WAM3 model, however, deviates from the dynamical model significantly. 

As in the artificial viscosity case, the angle-averaged wave action spectra as the func- 
tion of frequency exhibit distinct down-shift of the spectral maxima for both dynamical 
and Hasselmann equations (see Fig l32l) . Correspondence between time evolution of the 
amplitudes of the spectral maxima is also much better for WAM3 choice of damping, than 
for artificial viscosity case. 

Presumably, WAM3 damping term underestimate the effects of real damping at the very 
beginning of the evolution (when the effects of white capping are relatively important), 
and overestimates them on later stages of swell evolution. 
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5.3. Statistical experiments with WAM4 damping term 

The final third series of experiments have been done for the choice of WAM4 damping 
term. 

FigJ71-[T0] show temporal evolution of the total wave action, total energy, mean wave 
slope and mean wave frequency, which are divergent in time in this case. 

Figj33] show time evolution of angle-averaged wave action spectra as the functions of 
frequency for dynamical and Hasselmann equations. While as in the artificial viscosity 
and WAM3 cases we also observe distinct down-shift of the spectral maxima, the corre- 
spondence of the time evolution of the amplitudes of the spectral maxima is worse than 
in WAM3 case. 

This observation is especially surprising in view of the fact that historically WAM4 
damping has been invented as an improvement to WAM3 damping term. It is quite 
obvious that WAM4 damping is too strong for description of the reality at all stages of 
the swell evolution. 



6. Down-shift and angular spreading 

The major process of time-evolution of the swell is frequency down-shift. During T = 
933To the mean frequency has been decreased from cuq = 2 to cji = 0.6. On the last stage 
of the process, the mean frequency slowly decays as 

< >~ ^ (27) 

The Hasselmann equation has self-similar solution, describing the evolution of the swell 
n{k,t) (see [|13], [SS]). For this solution 

< >~ t-i/" (28) 

The difference between (^7^ and (1251) can be explained as follows. What we observed, is 
not a self-similar behavior. Indeed, a self-similarity presumes that the angular structure 
of the solution is constant in time. Meanwhile, we observed intensive angular spreading 
of the initially narrow in angle, almost one- dimensional wave spectrum. 

Level lines of the low-pass filtered dynamic and kinetic spectrum at the beginning of 
the simulation are presented on Fig J2Bll2Bl There are ten isolines on every figure. Values 
of level at maximum and minimum isolines are shown on every picture. Level lines of 
the low-pass filtered dynamic and kinetic spectrum for three different damping terms at 
the end of the simulation are presented on Fig j271l28l We observed development of 
bimodality in both experiments. This is in accordance with field observations [ [53] , [ [51] . 

One can see good correspondence between the results of both experiments. Comparison 
of time-evolution of the mean angular spreading, calculated from action and energy spectra 
is presented on Fig. [29l[30l 
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Figure 26. Level lines of the dynamic (left) and the kinetic (right) spectra Bit t — 67.1 Tq. 
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Dynamic equations 
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Figure 27. Level lines of the dynamic (left) and the kinetic (right) spectra ait — 3378.2 Tq. 
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Figure 28. (left)Level lines of the kinetic spectra in the WAM3 (left) WAM4 and damping 
cases at t = 3378.2 Tq. 
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Calculation from wave action spectrum 
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Figure 29. Average angular spreading (^J \9\n{k)dk^ / (^J n(k)dk^ as a function of time, 
calculated from wave action. Solid line - WAM3, dotted line - WAM4, dashed line - 
artificial viscosity, stars - dynamical equations 
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Calculation from wave energy spectrum 
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Figure 30. Average angular spreading ^/ \6\u!n{k)dkj / {j u!n{k)dkj as a function of time, 
calculated from wave energy. Solid line - WAM3, dotted line - WAM4, dashed hue - 
artificial viscosity, stars - dynamical equations 
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We see growing divergence between dynamic and kinetic models. But using WAM3 and 
WAM4 models leads to worse divergence. This is an additional argument against these 
variants of white-cap damping. 

One can expect that the angular spreading will be arrested at later times, and the 
spectra will take universal self-similar shape. 




Frequency 



Figure 31. Angle-averaged spectrum as a function of time for dynamical and Hasselmann 
equations for artificial viscosity case. 



7. Conclusion 

1. Our numerical experiment shows that the Hasselmann equation without any ad- 
ditional terms is applicable for description of gravity waves turbulence in the in- 
finite basin only if nonlinearity is low enough. The average wave steepness /i = 
{E^ /N'^)sqrt2E should be significantly less a certain critical level /xq = 0.095. To 
describe the wave turbulence of a moderate steepness 0.095 < /i < 0.15, one should 
add to the Hasselmann equation some additional term Sdiss modeling dissipation 
due to white-capping. So far there is no any theory making possible to find Sdiss 
analytically. This is a great challenge for hydro dynamicists. So far we can suggest 
that Sdiss depends on steepness very sharply. Hence the guess that white-capping is 
the threshold phenomenon formulated by Banner, Babanin, and Young [ [52] looks 
very plausible. We did not model too steep waves, but we can conjecture that if 
/i > 0.15, white capping is the leading nonlinear process and the weak-turbulent 
approach is not applicable at al. 
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Figure 32. Angle-averaged spectrum as a function of time for dynamical and Hasselmann 
equations a function of time for WAM3 case. 



2. Our results are not a surprise for designers of operational models for wave forecast- 
ing. They understood necessity to introduce into the models extra dissipative terms 
many years ago. Another story that the models of Sdiss routinely used in the WAM 
and WAVEWATCH models are too crude and do not grasp the most important 
feature of the white-capping — its threshold nature. In our opinion existing models 
of Sdiss overestimate dissipation at small values of steepness and underestimate it 
at large fi. Moreover, they overestimate dissipation in the area of the spectral peak 
and underestimate it in the spectral region of high wave numbers. We are planning 
to offer our own form of Sdiss extracted from our massive numerical simulation of 
Euler equation, but it will take some time and toil. 

3. We have to stress that we solved maximally idealized model. We did not take 
into account interaction of swell with atmosphere (see for instance [[55], [[56] ), 
interaction with ocean currents, deviation of the wave motion from potentiality as 
well as influence of turbulence in the surface-atmosphere boundary layer. In the 
future we plan to establish a contact with experimentalists to develop more realistic 
model of swell propagation in the real ocean. 

4. In one aspect our conclusions are very resolute. All existing experimental wave tanks 
cannot be used for modelling of the waves propagation in open sea. The mesoscopic 
effects in wave tanks are too strong for reasonable values of steepness. This pertains 
only to modeling of kinetics of energy containing region in the vicinity of spectral 
peak. The rear faces of spectral distributions, spectral tails, can be successfully 
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Figure 33. Angle-averaged spectrum as a function of time for dynamical and Hasselmann 
equations a function of time for WAM4 case. 



modelled. This was demonstrated by Toba in his classical work [ESj- Toba observed 
the Zakharov- Filonenko spectrum ~ uj~^ in the wave tank of moderate size. In 
our opinion this conclusion is very important. Some authors, trying to find the 
universal law of fetch dependence of mean energy and mean frequency put together 
data collected in ocean and in the experimental tanks. This mixed data hardly can 
be compared with any self sufficient theory. This question is discussed in details in 
the article [[57]. 

Another conclusion is more pessimistic. The results of numerical experiments show that it 
is very difficult to reproduce real ocean conditions in laboratory wave tank. Even the size 
of the tanks of 200 x 200 meters is not large enough to model ocean due to the presence 
of wave numbers grid discreteness. 
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A. "Forbes" list of 15 largest harmonics. 

Here one can find 15 largest harmonics at the end of calculations in the framework of 
dynamical equations. Average square of amplitudes in dissipation-less region was taken 
from smoothed spectrum, while all these harmonics exceed level |a^p = 1.4 x 10~^^. 
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B. Prom Dynamical Equations to Hasselmann Equation. 

Standard setup for numerical simulation of the dynamical equations (jlj) implies 27r x 27r 
domain in real space and gravity acceleration g = 1. Usage of the domain size equal 27r 
is convenient because in this case wave numbers are integers. 

In the contrary to dynamical equations, the kinetic equation f[T21) is formulated in 
terms of real physical variables and it is necessary to describe the transformation from 
the "dynamical" variables into to the "physical" ones. 

EqH] are invariant with respect to "stretching" transformation from "dynamical" to 
"real" variables: 

r]f = ar]!p,, k = —k', f = ar , g = vg\ (29) 



a 

t = J-t', L, = aL'^, Ly = aL'y (30) 



where prime denotes variables corresponding to dynamical equations. 

In current simulation we used the stretching coefficient a = 800.00, which allows 
to reformulate the statement of the problem in terms of real physics: we considered 
5026 m X 5026 m periodic boundary conditions domain of statistically uniform ocean with 
the same resolution in both directions and characteristic wave length of the initial condi- 
tion around 22 m. In oceanographic terms, this statement corresponds to the "duration- 
limited experiment". 
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